##############################
## Code to install packages ##
##############################
# install.packages(c(
#   "ggplot2",
#   "readxl",
#   "countrycode",
#   "ggrepel",
#   "ggpubr",
#   "pals",
#   "forcats",
#   "tidyverse",
#   "patchwork",
#   "conflicted",
#   "MASS"),
#   dep=T)




###################
## Load packages ##
###################
library(ggplot2)
library(readxl)
library(countrycode)
library(ggrepel)
library(ggpubr)
library(pals)
library(forcats)
library(tidyverse)
library(patchwork)
library(conflicted)
library(MASS)
conflict_prefer("select", "dplyr", "MASS")
conflict_prefer("filter", "dplyr", "stats")




###############
## Load data ##
###############
options(scipen=999)

## Set folder with replication files as working directory
setwd("C:/FOLDER")

dat <- read.csv("data.csv")




##############
## Figure 1 ##
##############
dt <- dat %>%
  filter(lc_elec_year==1 | uc_elec_year==1) %>%
  select(ctr_n, pspn_congruence, psns_congruence) %>%
  pivot_longer(!ctr_n) %>%
  mutate(name=factor(name, levels=c("pspn_congruence", "psns_congruence"), labels=c("Partisan", "Nationalization"))) %>%
  group_by(ctr_n) %>%
  mutate(ct=n()/2,
         ctr_n_ct=factor(paste0(ctr_n," (",ct,")")))

p <- ggplot(data=dt, aes(x=fct_rev(ctr_n_ct), y=value, fill=fct_rev(name))) +
  theme_bw() +
  theme(legend.position="bottom",
        legend.title=element_blank(),
        legend.text=element_text(size=12),
        axis.text.x=element_text(size=12, color="black"),
        axis.text.y=element_text(size=12, color="black"),
        axis.title.x=element_text(size=12),
        axis.title.y=element_blank(),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.caption=element_text(hjust=0, size=12),
        plot.caption.position="plot") +
  geom_boxplot(color="black", width=0.5, position=position_dodge(0.75), outlier.size=1) +
  ylab("Congruence") +
  scale_fill_manual(values=c("bisque3","royalblue4")) +
  guides(fill=guide_legend(nrow=1, ncol=3, reverse=T)) +
  coord_flip()

ggsave("Plots/Figure 1.png", width=6, height=6, dpi=300)




##############
## Figure 2 ##
##############
dt <- dat %>%
  filter(lc_elec_year==1 | uc_elec_year==1) %>%
  select(ctr_n, yr, pspn_congruence, psns_congruence) %>%
  rename(Partisan=pspn_congruence,
         Nationalization=psns_congruence) %>%
  mutate(avg=0)

dt_ct <- dt %>%
  group_by(ctr_n) %>%
  summarise(n=n()) %>%
  ungroup

dt_mean <- dt %>%
  group_by(ctr_n) %>%
  summarize_all(list(mean)) %>%
  mutate(avg=1)

dt2 <- dt %>%
  rbind(., dt_mean) %>%
  mutate(avg=factor(avg)) %>%
  mutate(cc=ifelse(avg==1, countrycode(ctr_n, "country.name", "iso3c"), "")) %>%
  mutate(cc=factor(cc)) %>%
  mutate(us_obs=factor(ifelse(ctr_n=="United States" & avg==0, 1, 0),
                       labels=c("Non-US", "US"))) %>%
  mutate(Nationalization=ifelse(cc=="ZWE", Nationalization+0.005, Nationalization),
         Nationalization=ifelse(cc=="MEX", Nationalization-0.005, Nationalization),
         Nationalization=ifelse(cc=="ARG", Nationalization-0.005, Nationalization),
         Nationalization=ifelse(cc=="NGA", Nationalization-0.007, Nationalization),
         Partisan=ifelse(cc=="AUS", Partisan-0.0015, Partisan),
         Partisan=ifelse(cc=="CHL", Partisan+0.0015, Partisan),
         Partisan=ifelse(cc=="ITA", Partisan-0.005, Partisan),
         Nationalization=ifelse(cc=="BRA", Nationalization-0.007, Nationalization))

p <- ggplot(data=subset(dt2, avg==0), aes(x=Partisan, y=Nationalization, shape=us_obs)) +
  theme_bw() +
  theme(legend.position="bottom",
        legend.title=element_blank(),
        legend.text=element_text(size=12),
        axis.text.x=element_text(size=12, color="black"),
        axis.text.y=element_text(size=12, color="black"),
        axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.caption=element_text(hjust=0, size=12),
        plot.caption.position="plot") +
  geom_point(aes(size=us_obs), alpha=0.2) +
  scale_shape_manual(values=c(19, 3)) +
  scale_size_manual(values=c(2,3)) +
  scale_x_continuous(breaks=seq(0,1,0.25), lim=c(0,1)) +
  scale_y_continuous(breaks=seq(0,1,0.25), lim=c(0,1)) +
  geom_text(data=subset(dt2, avg==1), aes(x=Partisan, y=Nationalization, label=cc), size=3.5)


ggsave("Plots/Figure 2.png", width=7, height=6, dpi=300)




##############
## Figure 3 ##
##############
ctry <- dat %>%
  select(ctr_n) %>%
  unique %>%
  mutate(id=seq(1,n(),1)) %>%
  mutate(ctr_n=str_replace_all(ctr_n, "Dominican Republic", "Dom. Rep."))

n_indvar <- 7
var_lab <- c("Symmetric",
             "DM Ratio",
             "Seat Ratio",
             "Time Diff",
             "Partial\nRenew",
             "Executive\nStrength",
             "Yrs Since\nPres")

cause_plot <- rbind(
  cbind(
    read.csv("Results/cause_plot_90.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      filter(`X..1`!="") %>%
      filter(`X.`!="_cons") %>%
      separate(X..1, into=c("lb_90", "ub_90"), sep=",", convert=TRUE),
    read.csv("Results/cause_plot_95.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      filter(`X..1`!="") %>%
      filter(`X.`!="_cons") %>%
      separate(X..1, into=c("lb_95", "ub_95"), sep=",", convert=TRUE) %>%
      select(lb_95, ub_95)
  ) %>%
    rename(var=1,
           pe=2) %>%
    mutate(model=c(rep(0, n_indvar), rep(1, n_indvar)),
           sample=0),
  cbind(
    read.csv("Results/cause_plot_nonus_90.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      filter(`X..1`!="") %>%
      filter(`X.`!="_cons") %>%
      separate(X..1, into=c("lb_90", "ub_90"), sep=",", convert=TRUE),
    read.csv("Results/cause_plot_nonus_95.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      filter(`X..1`!="") %>%
      filter(`X.`!="_cons") %>%
      separate(X..1, into=c("lb_95", "ub_95"), sep=",", convert=TRUE) %>%
      select(lb_95, ub_95)
  ) %>%
    rename(var=1,
           pe=2) %>%
    mutate(model=c(rep(0, n_indvar), rep(1, n_indvar)),
           sample=1)
) %>%
  mutate(pe=as.numeric(pe),
         sample=factor(sample,
                       labels=c("Full Sample", "Non-US Sample")),
         var=factor(var,
                    levels=unique(var),
                    labels=var_lab),
         model=factor(model,
                      levels=unique(model),
                      labels=c("Partisan", "Nationalization")),
         type=factor(paste0(model," (",sample,")"), levels=c("Partisan (Full Sample)",
                                                             "Partisan (Non-US Sample)",
                                                             "Nationalization (Full Sample)",
                                                             "Nationalization (Non-US Sample)")),
         sig=factor(ifelse(lb_90*ub_90>0, 1, 0)),
         pe=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", pe/10, pe),
         lb_90=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", lb_90/10, lb_90),
         ub_90=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", ub_90/10, ub_90),
         lb_95=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", lb_95/10, lb_95),
         ub_95=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", ub_95/10, ub_95))


p <- ggplot(cause_plot, aes(x=fct_rev(var), group=fct_rev(type), y=pe, color=type, shape=type, alpha=sig)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_text(size=10, color="black"),
        axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=10, color="black"),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=10),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.caption=element_text(hjust=0, size=10),
        plot.caption.position="plot",
        plot.margin=unit(c(5.5,12,5.5,10), "pt"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=12)) +
  geom_errorbar(aes(ymin=lb_90, ymax=ub_90, width=0), size=1, position=position_dodge(width=0.7), show.legend=F) +
  geom_errorbar(aes(ymin=lb_95, ymax=ub_95, width=0), size=0.5, position=position_dodge(width=0.7), show.legend=F) +
  geom_point(position=position_dodge(width=0.7), size=1.75) +
  geom_hline(yintercept=0, color="red") +
  scale_y_continuous(breaks=round(seq(-0.15,0.15,0.05),2),
                     lim=c(-0.15,0.15),
                     expand=c(0,0)) +
  scale_color_manual(name="legend",
                     values=c("royalblue4", "royalblue4", "bisque3", "bisque3")) +
  scale_shape_manual(name="legend",
                     values=c(19,17,19,17)) +
  scale_alpha_manual(values=c(0.35,1), guide="none") +
  facet_wrap(~model, ncol=1, scales="free_y", strip.position="top") +
  guides(color=guide_legend(nrow=2, ncol=2)) +
  coord_flip()

ggsave(paste0("Plots/Figure 3.png"), width=8, height=8, dpi=300)




##############
## Figure 4 ##
##############
## Figure 4A
consq_plot <- cbind(
  read.csv("Results/consq_90.csv") %>%
    mutate_all(funs(str_replace(., "=", ""))) %>%
    rename(var=1,
           pe=2) %>%
    filter(`X..1`!="") %>%
    separate_wider_delim(`X..1`, ",", names=c("lb_90", "ub_90")),
  read.csv("Results/consq_95.csv") %>%
    mutate_all(funs(str_replace(., "=", ""))) %>%
    rename(var=1,
           pe=2) %>%
    filter(`X..1`!="") %>%
    separate_wider_delim(`X..1`, ",", names=c("lb_95", "ub_95")) %>%
    select(contains("95"))
) %>%
  mutate_at(vars(pe, lb_90, ub_90, lb_95, ub_95), as.numeric) %>%
  mutate(var=factor(var,
                    levels=unique(var),
                    labels=c("Not Symmetric",
                             "Symmetric"))) %>%
  mutate(sig=factor(ifelse((lb_90*ub_90)<0,0,1), levels=c(0,1))) %>%
  mutate(ind_var=factor(c(rep("Partisan", 2), rep("Nationalization", 2))))

p4a <- ggplot(consq_plot, aes(x=var, y=pe, color=fct_rev(ind_var), shape=fct_rev(ind_var), alpha=sig)) +
  theme(axis.title.y=element_text(size=12, color="black"),
        axis.text.y=element_text(size=12, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.title=element_text(size=14, hjust=0.5),
        plot.caption=element_text(hjust=0, size=10),
        plot.caption.position="plot",
        plot.margin=unit(c(5,15,5,5), "pt")) +
  ggtitle("(A) Average Marginal Effects") +
  geom_errorbar(aes(ymin=lb_90, ymax=ub_90, width=0), size=1, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(ymin=lb_95, ymax=ub_95, width=0), size=0.5, position=position_dodge(width=0.4)) +
  geom_point(size=1.75, position=position_dodge(width=0.4)) +
  geom_hline(yintercept=0, linetype="solid", color="red") +
  scale_color_manual(values=c("royalblue4","bisque3"), drop=F) +
  scale_shape_manual(values=c(19, 17)) +
  scale_alpha_manual(values=c(1,1), guide="none") +
  scale_y_continuous(breaks=seq(-2, 2, 1),
                     lim=c(-2,2),
                     name="AME",
                     sec.axis=sec_axis(~ (.+2)*100, name="Number of Observations")) +
  theme(strip.background=element_blank(),
        strip.text=element_text(size=10),
        strip.placement="outside",
        panel.spacing=unit(0.25, "cm")) +
  geom_segment(aes(x=1, xend=1, y=-2, yend=-0.73),
               linewidth=10,
               color="grey",
               alpha=0.2,
               position=position_nudge(x=0.00)) +
geom_segment(aes(x=2, xend=2, y=-2, yend=-0.16),
             linewidth=10,
             color="grey",
             alpha=0.2,
             position=position_nudge(x=0.00))



## Figure 4B
response_curve <- function(t, lag.coef, x.coef){
  response <- list(temporary.shock=rep(0,length(t)),
                   permanent.shock=rep(0,length(t)))
  response$temporary.shock <- x.coef * lag.coef^(t-1)
  response$permanent.shock <- cumsum(response$temporary.shock)
  return(response)
}

lr_sim <- function(T, S, pspn){

  coef <- read.csv("Results/consq_coef.csv") %>%
    mutate_all(funs(str_replace(., "=", ""))) %>%
    mutate_all(as.numeric) %>%
    t() %>%
    as.data.frame %>%
    filter(!is.na(V1) & V1!=0) %>%
    pull(V1)

  vcov <- read.csv("Results/consq_vcov.csv") %>%
    mutate_all(funs(str_replace(., "=", ""))) %>%
    mutate_all(as.numeric) %>%
    select(-1) %>%
    filter(rowSums(.)!=0) %>%
    select_if(~sum(.)!=0) %>%
    as.matrix


  t <- 1:T
  sim.coef <- mvrnorm(S, coef, vcov)
  sim.temp <- matrix(NA, S, length(t))

  if (pspn==1){

    response <- response_curve(t, coef[1], (coef[3] + coef[4])*0.1373)

    for(i in 1:S){
      sim.lag <- sim.coef[i,1]
      sim.x <- (sim.coef[i,3] + sim.coef[i,4])*0.1373
      sim.shock <- response_curve(t, sim.lag, sim.x)
      sim.temp[i,] <- sim.shock$permanent.shock
    }

  } else{

    response <- response_curve(t, coef[1], (coef[5] + coef[6])*0.1917)

    for(i in 1:S){
      sim.lag <- sim.coef[i,1]
      sim.x <- (sim.coef[i,5] + sim.coef[i,6])*0.1917
      sim.shock <- response_curve(t, sim.lag, sim.x)
      sim.temp[i,] <- sim.shock$permanent.shock
    }

  }

  se_sim <- apply(sim.temp, 2, sd)

  out <- data.frame(pe=response$permanent.shock,
                    lb_90=response$permanent.shock - se_sim*1.645,
                    ub_90=response$permanent.shock + se_sim*1.645,
                    lb_95=response$permanent.shock - se_sim*1.96,
                    ub_95=response$permanent.shock + se_sim*1.96)

  return(out)

}


set.seed(2024)
lr_plot <- rbind(
  lr_sim(5, 1000, 1),
  lr_sim(5, 1000, 2)
) %>%
  mutate(type=factor(c(rep(0, 5), rep(1, 5)),
                     labels=c("Partisan", "Nationalization"))) %>%
  mutate(t=rep(seq(1,5,1), 2))

p4b <- ggplot(lr_plot, aes(x=t, y=pe, color=type, shape=type)) +
  theme(axis.title.y=element_text(size=12, color="black"),
        axis.text.y=element_text(size=12, color="black"),
        axis.title.x=element_text(size=12, color="black"),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=10),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.title=element_text(size=14, hjust=0.5),
        plot.caption=element_text(hjust=0, size=12),
        plot.caption.position="plot",
        plot.margin=unit(c(5,5,5,15), "pt")) +
  ggtitle("(B) Short-Run Response") +
  geom_errorbar(aes(ymin=lb_90, ymax=ub_90, width=0), size=1) +
  geom_errorbar(aes(ymin=lb_95, ymax=ub_95, width=0), size=0.5) +
  geom_point(size=1.75) +
  geom_hline(yintercept=0, linetype="solid", color="red") +
  scale_color_manual(values=c("royalblue4","bisque3"), drop=F) +
  scale_shape_manual(values=c(19, 17)) +
  xlab("Years after One Standard Deviation Increase in Congruence") +
  ylab("Estimated Change in Gov't Spending")


## Combine and save
p4 <- ggarrange(p4a, p4b, nrow=1, common.legend=T, legend="bottom")

ggsave(p4, filename=paste0("Plots/Figure 4.png"), width=12, height=4, dpi=300, background="white")




###############
## Figure A1 ##
###############
pd <- dat %>%
  filter(lc_elec_year==1 | uc_elec_year==1) %>%
  select(ctr_n, yr, pspn_congruence, psns_congruence) %>%
  arrange(ctr_n) %>%
  mutate(ctr_n=factor(ctr_n, levels=unique(ctr_n))) %>%
  gather(key="var", value="value", -c(ctr_n, yr)) %>%
  mutate(var=factor(var, levels=c("pspn_congruence", "psns_congruence"), labels=c("Partisan", "Nationalization")))

p <- ggplot(pd, aes(x=yr, y=value)) +
  xlab("Year") +
  ylab("Congruence") +
  scale_x_continuous(limits=c(min(pd$yr)-1,max(pd$yr)+1)) +
  ylim(0.1,1) +
  theme_bw() +
  theme(legend.position="bottom",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=10),
        plot.title=element_text(size=14, hjust=0.5),
        axis.text.x=element_text(size=10, color="black"),
        axis.text.y=element_text(size=10, color="black"),
        axis.title.x=element_text(size=10),
        axis.title.y=element_text(size=10),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=12)) +
  geom_line(aes(color=var), size=0.5) +
  geom_point(aes(color=var, shape=var), size=1.25) +
  scale_color_manual(values=c("royalblue4", "bisque3")) +
  scale_x_continuous(breaks=seq(1920,2020,25)) +
  facet_wrap(~ctr_n, ncol=3) +
  theme(strip.text=element_text(size=10))

ggsave("Plots/Figure A1.png", width=8, height=10)




###############
## Figure C1 ##
###############
consq_plot <- rbind(
  cbind(
    read.csv("Results/consq_90.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      rename(var=1,
             pe=2) %>%
      filter(`X..1`!="") %>%
      separate_wider_delim(`X..1`, ",", names=c("lb_90", "ub_90")),
    read.csv("Results/consq_95.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      rename(var=1,
             pe=2) %>%
      filter(`X..1`!="") %>%
      separate_wider_delim(`X..1`, ",", names=c("lb_95", "ub_95")) %>%
      select(contains("95"))
  ),
  cbind(
    read.csv("Results/consq_cbi_90.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      rename(var=1,
             pe=2) %>%
      filter(`X..1`!="") %>%
      separate_wider_delim(`X..1`, ",", names=c("lb_90", "ub_90")),
    read.csv("Results/consq_cbi_95.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      rename(var=1,
             pe=2) %>%
      filter(`X..1`!="") %>%
      separate_wider_delim(`X..1`, ",", names=c("lb_95", "ub_95")) %>%
      select(contains("95"))
  ),
  cbind(
    read.csv("Results/consq_leftgov_90.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      rename(var=1,
             pe=2) %>%
      filter(`X..1`!="") %>%
      separate_wider_delim(`X..1`, ",", names=c("lb_90", "ub_90")),
    read.csv("Results/consq_leftgov_95.csv") %>%
      mutate_all(funs(str_replace(., "=", ""))) %>%
      rename(var=1,
             pe=2) %>%
      filter(`X..1`!="") %>%
      separate_wider_delim(`X..1`, ",", names=c("lb_95", "ub_95")) %>%
      select(contains("95"))
  )) %>%
  mutate_at(vars(pe, lb_90, ub_90, lb_95, ub_95), as.numeric) %>%
  mutate(var=factor(var,
                    levels=unique(var),
                    labels=c("Not Symmetric",
                             "Symmetric"))) %>%
  mutate(sig=factor(ifelse((lb_90*ub_90)<0,0,1), levels=c(0,1))) %>%
  mutate(ind_var=factor(rep(c(rep("Partisan", 2), rep("Nationalization", 2)), 3))) %>%
  mutate(model=factor(sort(rep(seq(1,3,1), 4)),
                      labels=c("Model 1", "Model 2", "Model 3")))



p <- ggplot(consq_plot, aes(x=var, y=pe, color=fct_rev(ind_var), shape=model)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_text(size=8, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=10, color="black"),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=10),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.caption=element_text(hjust=0, size=10),
        plot.caption.position="plot") +
  geom_errorbar(aes(ymin=lb_90, ymax=ub_90, width=0), size=1, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(ymin=lb_95, ymax=ub_95, width=0), size=0.5, position=position_dodge(width=0.4)) +
  geom_point(size=1.75, position=position_dodge(width=0.4)) +
  geom_hline(yintercept=0, linetype="solid", color="red") +
  scale_color_manual(values=c("royalblue4","bisque3"), drop=F) +
  scale_alpha_manual(values=c(0.3,1), guide="none") +
  theme(strip.background=element_blank(),
        strip.text=element_text(size=10),
        strip.placement="outside",
        panel.spacing=unit(0.25, "cm")) +
  labs(caption="Note: thick and thin bars show 90% and 95% confindence intervals, respectively. Full results\ncan be found in the appendix.")

ggsave(paste0("Plots/Figure C1.png"), width=6, height=4)




###############
## Figure D1 ##
###############
ct_list <- sort(ctry$ctr_n)

n_obs <- read.csv("Results/cause_jk_plot_90.csv") %>%
  mutate_all(funs(str_replace(., "=", ""))) %>%
  filter(`X.`=="N") %>%
  select(-`X.`) %>%
  unlist() %>%
  as.numeric() %>%
  na.omit()

cause_plot <- cbind(
  read.csv("Results/cause_jk_plot_90.csv") %>%
    mutate_all(funs(str_replace(., "=", ""))) %>%
  filter(`X..1`!="") %>%
  filter(`X.`!="_cons") %>%
  rename_with(~c("var",
                 unlist(mapply(function(pe, se) c(pe, se),
                               paste0("pe", 1:17),
                               paste0("ci", 1:17))))) %>%
  pivot_longer(
    cols=starts_with("pe") | starts_with("ci"),
    names_to=c(".value", "ct"),
    names_pattern="(pe|ci)(\\d+)"
  ) %>%
  mutate(ct=as.numeric(ct)) %>%
  arrange(ct) %>%
  separate(ci, into=c("lb_90", "ub_90"), sep=",", convert=TRUE),
  read.csv("Results/cause_jk_plot_95.csv") %>%
    mutate_all(funs(str_replace(., "=", ""))) %>%
    filter(`X..1`!="") %>%
    filter(`X.`!="_cons") %>%
    rename_with(~c("var",
                   unlist(mapply(function(pe, se) c(pe, se),
                                 paste0("pe", 1:17),
                                 paste0("ci", 1:17))))) %>%
    pivot_longer(
      cols=starts_with("pe") | starts_with("ci"),
      names_to=c(".value", "ct"),
      names_pattern="(pe|ci)(\\d+)"
    ) %>%
    mutate(ct=as.numeric(ct)) %>%
    arrange(ct) %>%
    separate(ci, into=c("lb_95", "ub_95"), sep=",", convert=TRUE) %>%
    select(lb_95, ub_95)
) %>%
  mutate(ct=factor(ct, labels=paste0(ct_list, " (N=", n_obs, ")")),
         pe=as.numeric(pe),
         model=rep(c(rep(0, n_indvar), rep(1, n_indvar)), length(ct_list)),
         var=factor(var,
                    levels=unique(var),
                    labels=var_lab),
         model=factor(model,
                      levels=unique(model),
                      labels=c("Partisan", "Nationalization")),
         pe=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", pe/10, pe),
         lb_90=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", lb_90/10, lb_90),
         ub_90=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", ub_90/10, ub_90),
         lb_95=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", lb_95/10, lb_95),
         ub_95=ifelse(var=="Partial\nRenew" | var=="Executive\nStrength", ub_95/10, ub_95))


## Figure D1a
p <- ggplot(cause_plot %>%
           filter(model=="Partisan"), aes(x=fct_rev(var), group=fct_rev(ct), y=pe, color=ct)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_text(size=12, color="black"),
        axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="right",
        legend.title=element_text(size=12),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.caption=element_text(hjust=0, size=12),
        plot.caption.position="plot",
        plot.margin=unit(c(5.5,12,5.5,10), "pt"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=12)) +
  geom_errorbar(aes(ymin=lb_90, ymax=ub_90, width=0), size=1, position=position_dodge(width=0.9), show.legend=F) +
  geom_errorbar(aes(ymin=lb_95, ymax=ub_95, width=0), size=0.5, position=position_dodge(width=0.9), show.legend=F) +
  geom_point(position=position_dodge(width=0.9), size=1.25, shape=16) +
  geom_hline(yintercept=0, color="red") +
  scale_y_continuous(breaks=round(seq(-0.15,0.15,0.05),2),
                     lim=c(-0.15,0.175),
                     expand=c(0,0)) +
  scale_color_manual(values=setNames(pals::polychrome(n=length(ct_list)), paste0(ct_list, " (N=", n_obs, ")"))) +
    facet_wrap(~model, ncol=1, scales="free_y", strip.position="top") +
  guides(color=guide_legend(ncol=1)) +
  coord_flip() +
  labs(color="Excluded Country:",
       caption="Note: thick and thin bars show 90% and 95% confindence intervals, respectively. The numbers in the legend specify the number of\nobservations in the sample after a given country is excluded. Estimates for the Partial Renew and Executive Strength variables have\nbeen rescaled for this plot by dividing by 10.")

ggsave(paste0("Plots/Figure D1a.png"), width=10, height=10)



## Figure D1b
p <- ggplot(cause_plot %>%
              filter(model=="Nationalization"), aes(x=fct_rev(var), group=fct_rev(ct), y=pe, color=ct)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_text(size=12, color="black"),
        axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="right",
        legend.title=element_text(size=12),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.caption=element_text(hjust=0, size=12),
        plot.caption.position="plot",
        plot.margin=unit(c(5.5,12,5.5,10), "pt"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=12)) +
  geom_errorbar(aes(ymin=lb_90, ymax=ub_90, width=0), size=1, position=position_dodge(width=0.9), show.legend=F) +
  geom_errorbar(aes(ymin=lb_95, ymax=ub_95, width=0), size=0.5, position=position_dodge(width=0.9), show.legend=F) +
  geom_point(position=position_dodge(width=0.9), size=1.25, shape=16) +
  geom_hline(yintercept=0, color="red") +
  scale_y_continuous(breaks=round(seq(-0.15,0.15,0.05),2),
                     lim=c(-0.15,0.175),
                     expand=c(0,0)) +
  scale_color_manual(values=setNames(pals::polychrome(n=length(ct_list)), paste0(ct_list, " (N=", n_obs, ")"))) +
  facet_wrap(~model, ncol=1, scales="free_y", strip.position="top") +
  guides(color=guide_legend(ncol=1)) +
  coord_flip() +
  labs(color="Excluded Country:",
       caption="Note: thick and thin bars show 90% and 95% confindence intervals, respectively. The numbers in the legend specify the number of\nobservations in the sample after a given country is excluded. Estimates for the Partial Renew and Executive Strength variables have\nbeen rescaled for this plot by dividing by 10.")

ggsave(paste0("Plots/Figure D1b.png"), width=10, height=10)




###############
## Figure D2 ##
###############
ct_list2 <- ct_list[!ct_list %in% c("Liberia", "Zimbabwe")]

n_obs2 <- read.csv("Results/consq_jk_90.csv") %>%
  mutate_all(funs(str_replace(., "=", ""))) %>%
  filter(`X.`=="N") %>%
  select(-`X.`) %>%
  select(where(~ !any(.=="311"))) %>%
  unlist() %>%
  as.numeric() %>%
  na.omit()

consq_plot <- cbind(
  read.csv("Results/consq_jk_90.csv") %>%
    mutate_all(funs(str_replace(., "=", ""))) %>%
    filter(`X..1`!="") %>%
    rename_with(~c("var",
                   unlist(mapply(function(pe, se) c(pe, se),
                                 paste0("pe", 1:17),
                                 paste0("ci", 1:17))))) %>%
    mutate(var=c("pspn_0", "pspn_1", "psns_0", "psns_1")) %>%
    pivot_longer(
      cols=starts_with("pe") | starts_with("ci"),
      names_to=c(".value", "ct"),
      names_pattern="(pe|ci)(\\d+)"
    ) %>%
    separate(ci, into=c("lb_90", "ub_90"), sep=",", convert=TRUE),
  read.csv("Results/consq_jk_95.csv") %>%
    mutate_all(funs(str_replace(., "=", ""))) %>%
    filter(`X..1`!="") %>%
    rename_with(~c("var",
                   unlist(mapply(function(pe, se) c(pe, se),
                                 paste0("pe", 1:17),
                                 paste0("ci", 1:17))))) %>%
    mutate(var=c("pspn_0", "pspn_1", "psns_0", "psns_1")) %>%
    pivot_longer(
      cols=starts_with("pe") | starts_with("ci"),
      names_to=c(".value", "ct"),
      names_pattern="(pe|ci)(\\d+)"
    ) %>%
    separate(ci, into=c("lb_95", "ub_95"), sep=",", convert=TRUE) %>%
    select(lb_95, ub_95)
) %>%
  mutate(ct=as.numeric(ct)) %>%
  filter(ct!=10 & ct!=17) %>%
  mutate(ct=factor(ct, labels=paste0(ct_list2, " (N=", n_obs2, ")"))) %>%
  arrange(ct) %>%
  mutate_at(vars(pe, lb_90, ub_90, lb_95, ub_95), as.numeric) %>%
  mutate(type=factor(case_when(
    grepl("_1", var) ~ "Symmetric",
    grepl("_0", var) ~ "Not Symmetric"
  )),
  var=factor(case_when(
    grepl("pspn", var) ~ "Partisan",
    grepl("psns", var) ~ "Nationalization"
  )))


## Figure D2a
p <- ggplot(consq_plot %>%
           filter(var=="Partisan"), aes(x=type, y=pe, group=ct, color=ct)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_text(size=12, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="right",
        legend.title=element_text(size=12, color="black"),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.title=element_text(size=14, hjust=0.5),
        plot.caption=element_text(hjust=0, size=10),
        plot.caption.position="plot") +
  geom_errorbar(aes(ymin=lb_90, ymax=ub_90, width=0), size=1, position=position_dodge(width=0.6)) +
  geom_errorbar(aes(ymin=lb_95, ymax=ub_95, width=0), size=0.5, position=position_dodge(width=0.6)) +
  geom_point(size=1.75, position=position_dodge(width=0.6)) +
  geom_hline(yintercept=0, linetype="solid", color="red") +
  scale_color_manual(values=c("royalblue4","bisque3"), drop=F) +
  scale_shape_manual(values=c(19, 17)) +
  scale_alpha_manual(values=c(1,1), guide="none") +
  scale_y_continuous(breaks=seq(-2, 2, 1),
                     lim=c(-2.25,2.5),
                     name="AME") +
  theme(strip.background=element_blank(),
        strip.text=element_text(size=10),
        strip.placement="outside",
        panel.spacing=unit(0.25, "cm")) +
    scale_color_manual(values=setNames(pals::polychrome(n=length(ct_list2)), paste0(ct_list2, " (N=", n_obs2, ")"))) +
    labs(color="Excluded Country:",
         caption="Note: thick and thin bars show 90% and 95% confindence intervals, respectively. The numbers in the legend specify the number\nof observations in the sample after a given country is excluded.")

ggsave(paste0("Plots/Figure D2a.png"), width=8, height=4)


## Figure D2b
p <- ggplot(consq_plot %>%
              filter(var=="Nationalization"), aes(x=type, y=pe, group=ct, color=ct)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_text(size=12, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="right",
        legend.title=element_text(size=12, color="black"),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey"),
        plot.title=element_text(size=14, hjust=0.5),
        plot.caption=element_text(hjust=0, size=10),
        plot.caption.position="plot") +
  geom_errorbar(aes(ymin=lb_90, ymax=ub_90, width=0), size=1, position=position_dodge(width=0.6)) +
  geom_errorbar(aes(ymin=lb_95, ymax=ub_95, width=0), size=0.5, position=position_dodge(width=0.6)) +
  geom_point(size=1.75, position=position_dodge(width=0.6)) +
  geom_hline(yintercept=0, linetype="solid", color="red") +
  scale_color_manual(values=c("royalblue4","bisque3"), drop=F) +
  scale_shape_manual(values=c(19, 17)) +
  scale_alpha_manual(values=c(1,1), guide="none") +
  scale_y_continuous(breaks=seq(-2, 2, 1),
                     lim=c(-2.25,2.5),
                     name="AME") +
  theme(strip.background=element_blank(),
        strip.text=element_text(size=10),
        strip.placement="outside",
        panel.spacing=unit(0.25, "cm")) +
  scale_color_manual(values=setNames(pals::polychrome(n=length(ct_list2)), paste0(ct_list2, " (N=", n_obs2, ")"))) +
  labs(color="Excluded Country:",
       caption="Note: thick and thin bars show 90% and 95% confindence intervals, respectively. The numbers in the legend specify the number\nof observations in the sample after a given country is excluded.")

ggsave(paste0("Plots/Figure D2b.png"), width=8, height=4)